#include "plant_1DSOS_RDM_v2.5.h"
// --------------------------------------------------
plant_1DSOS_RDM::plant_1DSOS_RDM(double idt,int iDim,int iL,int iMode, double iE_D,double iE_S,double iE_N)
{
	// Options
	Dim = iDim;
	if(Dim == 1){
		Lx = iL;
		Ly = 1;
		N_CL = 3;
	}else{
		Lx = iL;
		Ly = iL;
		N_CL = 5;
	}

	mode = iMode;
	assert(2*iMode<=iL);
	E_D = iE_D;
	E_S = iE_S;
	E_N = iE_N;

	// 2 Inputs: #0 = T, #1 = W
	int input_num = 2;
	int input_dim[2] = {1,1};
	// 4 outputs:
	//  Output #0: h
	//  Output #1: r2
	//  Output #2: m2
	//  Output #3: a
	int output_num = 4;
	int output_dim[4];
	output_dim[0] = 1;
	output_dim[1] = 1;
	output_dim[2] = 1;
	output_dim[3] = 1;
	block::initialize(input_num,input_dim,output_num,output_dim,idt);
	
	// States
	lattice = new int*[Ly];
	lattice_stat =  new int*[Ly];
	xcov_val = new double[Lx];
	for(int i=0;i<Ly;i++){
		lattice[i] = new int[Lx];
		lattice_stat[i] = new int[Lx];
	}
	p_migr = new double[Ly];
	p_desorp = new double[Ly];
	// Initialize random number generator
	gsl_rng_env_setup();
	gsl_T = gsl_rng_default;
	gsl_r = gsl_rng_alloc(gsl_T);
	gsl_rng_set(gsl_r,0);

	double SQRT_TWO = sqrt(2.0);
	double SQRT_PI = sqrt(pi);
	reset();
	ofstream flog("plant_1DSOS_RDM_v2.log");
	flog << "Dim     = " << Dim << "\n";
	flog << "Lx      = " << Lx << "\n";
	flog << "Ly      = " << Ly << "\n";
	flog << "E_D     = " << E_D << "\n";
	flog << "E_S     = " << E_S << "\n";
	flog << "E_N     = " << E_N << "\n";
	flog.close();
}
// --------------------------------------------------
plant_1DSOS_RDM::~plant_1DSOS_RDM()
{
	int i;
	for(i=0;i<Ly;i++)
	{
		delete []lattice[i];
		delete []lattice_stat[i];
	}
	delete []lattice;
	delete []lattice_stat;
	delete []xcov_val;
	delete []p_migr;
	delete []p_desorp;
	gsl_rng_free(gsl_r);
}
// --------------------------------------------------
void plant_1DSOS_RDM::reset()
{
	block::time = 0.0;

	// Inputs
	T = 0;
	W = 0;

	// States
	int x,y;
	for(y=0;y<Ly;y++){
		for(x=0;x<Lx;x++){
			lattice[y][x] = 0;
		}
	}
	for(y=0;y<Ly;y++){
		for(x=0;x<Lx;x++){
			init_local_stat(x,y);
		}
	}

	// Outputs
	h  = 0.0;
	r2 = 0.0;
	m2 = 0.0;
	a  = Lx;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::init_local_stat(int x,int y){
	if(Dim == 2){
		lattice_stat[y][x]=(lattice[y][(x+Lx-1)%Lx]>=lattice[y][x])
		+(lattice[y][(x+Lx+1)%Lx]>=lattice[y][x])
		+(lattice[(y+Ly-1)%Ly][x]>=lattice[y][x])
		+(lattice[(y+Ly+1)%Ly][x]>=lattice[y][x]);
	}
	else{
		lattice_stat[y][x]=(lattice[y][(x+Lx-1)%Lx]>=lattice[y][x])+(lattice[y][(x+Lx+1)%Lx]>=lattice[y][x]);
	}
	return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::init_migr_rate()
{
	int i;
	double k_m;

	k_m=2*T*KB2H;

	for(i=0;i<4;i++)
	{
		MIGR_RT[i]=k_m*exp(-(E_S+i*E_N)*EV2KB/(T+0.0));
	}
    MIGR_RT[4]=0;
	return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::init_desorp_rate(){
	int i;
	double k_m;

	k_m=2*T*KB2H;

	for(i=0;i<5;i++)
	{
		//DESORP_RT[i]=k_m*exp(-(E_D+i*E_N)*EV2KB/(temperature+0.0));
		DESORP_RT[i]=0;
	}
	return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::init_p_migr()
{
	int x,y;

	migr_ttl=0.0;
	for(y=0;y<Ly;y++)
	{
		p_migr[y]=0.0;
		for(x=0;x<Lx;x++)
		{
			p_migr[y]+=MIGR_RT[lattice_stat[y][x]];
			
		}
		migr_ttl+=p_migr[y];
	}
	return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::init_p_desorp()
{
	int x,y;

	desorp_ttl=0.0;
	for(y=0;y<Ly;y++)
	{
		p_desorp[y]=0.0;
		for(x=0;x<Lx;x++)
		{
			p_desorp[y]+=DESORP_RT[lattice_stat[y][x]];
			
		}
		desorp_ttl+=p_desorp[y];
	}
	return;
}
// ---------------------------------------------------
void plant_1DSOS_RDM::update(double sysTime,double iT,double iW)
{
	if(time>sysTime)
		return;
	
	double rdm;
	while(time<=sysTime)
	{
		if(iT!=T)
			set_temperature(iT);
		if(iW!=W)
			set_ads_rate(iW);
		
		rate_ttl = ads_ttl+migr_ttl+desorp_ttl;
		// rdm = RandRand()*rate_ttl;
		rdm = gsl_rng_uniform(gsl_r)*rate_ttl;
		if(rdm<ads_ttl)
			adsorption();
		else if(rdm<=migr_ttl+ads_ttl)
			migration();
		else
			desorption();

		time += OneMCTime(rate_ttl);
	}
	outputs[0][0] = getH();
	outputs[1][0] = getR2();
	outputs[2][0] = getM2();
	outputs[3][0] = getAutoCorLength();
	return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::update(double sysTime)
{
	update(sysTime,inputs[0][0],inputs[1][0]);
}
// --------------------------------------------------
void plant_1DSOS_RDM::getSnapshot(const string FileName)
{
	ofstream fout(FileName.c_str());
	fout << "Lx   = " << Lx << "\n";
	fout << "Ly   = " << Ly << "\n";
	fout << "Time = " << time << "\n";
	int x,y;
	for(y=0;y<Ly;y++){
		for(x=0;x<Lx;x++){
			fout << setw(5) << lattice[y][x];
			if(x!=Lx-1){fout << "\n";}
		}
		if(y!=Ly-1){fout << "\n";}
	}
	fout.close();
}
// --------------------------------------------------
double plant_1DSOS_RDM::getH()
{
	int x,y;
	double h_avg,sum = 0;
	for(y=0;y<Ly;y++)
	{
		for(x=0;x<Lx;x++)
		{
			sum+=lattice[y][x];
		}
	}
	h_avg = sum/Ly/Lx;
	return h_avg;
}
// --------------------------------------------------
double plant_1DSOS_RDM::getM2(){
	int x,y;
	double sum = 0;
	double m2;
	for(y=0;y<Ly;y++)
	{
		for(x=0;x<Lx-1;x++)
		{
			sum+=(lattice[y][x]-lattice[y][x+1])*(lattice[y][x]-lattice[y][x+1]);
		}
		sum+=(lattice[y][Lx-1]-lattice[y][0])*(lattice[y][Lx-1]-lattice[y][0]);
	}
	m2 = sum/Ly/Lx;
	return m2;
}
// --------------------------------------------------
double plant_1DSOS_RDM::getR2(){
	int x,y;
	double R2;
	double sum=0;
	double h_avg;
	
	h_avg = getH();
	for(y=0;y<Ly;y++){
		for(x=0;x<Lx;x++)		{
			sum+=(lattice[y][x]-h_avg)*(lattice[y][x]-h_avg);
		}
	}
	R2 = sum/Ly/Lx;
	return R2;
}
// --------------------------------------------------
double plant_1DSOS_RDM::getAutoCorLength(){
	int i,m;
	double h_avg = getH();
	double xcov_val_max=0.0;
	m = 0;
	for(i=0;i<Lx-1-m;i++){
		xcov_val_max += (lattice[0][i]-h_avg)*(lattice[0][i+m]-h_avg);
 	}
	
	double xcov_val;
	for(m=1;m<Lx-1;m++){
		xcov_val = 0.0;
		for(i=0;i<Lx-1-m;i++){
			xcov_val+=(lattice[0][i]-h_avg)*(lattice[0][i+m]-h_avg);
		}
		if(xcov_val<0.3679*xcov_val_max){
			break;
		}
	}
	return m*0.2e-9;
}
// --------------------------------------------------
void plant_1DSOS_RDM::set_temperature(double T_set)
{
	T = T_set;
	init_migr_rate();
	init_desorp_rate();
	init_p_migr();
	init_p_desorp();
	k_m2d = exp(-(E_D-E_S)*EV2KB/(T+0.0));
	return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::set_ads_rate(double W_set)
{
	W=W_set;
	ads_ttl=W*Lx*Ly;
	return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::add_remove_atom(int x, int y, int num){
	int xpos[5],ypos[5];
	double site_rate_mig, site_rate_desorp;
	int i;

	GetNeighbors(x,y,xpos,ypos);
	for(i=0;i<N_CL;i++){
		site_rate_mig = MIGR_RT[lattice_stat[ypos[i]][xpos[i]]];
		site_rate_desorp = DESORP_RT[lattice_stat[ypos[i]][xpos[i]]];
		migr_ttl  -= site_rate_mig;
		desorp_ttl -= site_rate_desorp;
		p_migr[ypos[i]] -= site_rate_mig;
		p_desorp[ypos[i]] -= site_rate_desorp;
	}

	if((lattice[y][x]==0) && (num==-1))
	{}
	else{
		lattice[y][x]+=num;
	}

	for(i=0;i<N_CL;i++){
		init_local_stat(xpos[i],ypos[i]);
	}

	for(i=0;i<N_CL;i++){
		site_rate_mig=MIGR_RT[lattice_stat[ypos[i]][xpos[i]]];
		site_rate_desorp=DESORP_RT[lattice_stat[ypos[i]][xpos[i]]];
		migr_ttl+=site_rate_mig;
		desorp_ttl+=site_rate_desorp;
		p_migr[ypos[i]]+=site_rate_mig;
		p_desorp[ypos[i]]+=site_rate_desorp;
	}
	return;
}
// --------------------------------------------------	
void plant_1DSOS_RDM::GetNeighbors(int x, int y, int *xpos, int *ypos)
{
	if(Dim==2){
		xpos[0]=x;
		xpos[1]=x-1;
		xpos[2]=x;
		xpos[3]=x+1;
		xpos[4]=x;
		if(xpos[1]<0) xpos[1]+=Lx;
		if(xpos[3]>=Lx) xpos[3]-=Lx;

		ypos[0]=y;
		ypos[1]=y;
		ypos[2]=y-1;
		ypos[3]=y;
		ypos[4]=y+1;
		if(ypos[2]<0) ypos[2]+=Ly;
		if(ypos[4]>=Ly) ypos[4]-=Ly;
	}
	else{
		xpos[0] = x;
		xpos[1] = x-1;
		xpos[2] = x+1;
		if(xpos[1]<0) xpos[1]+=Lx;
		if(xpos[2]>=Lx) xpos[2]-=Lx;
		ypos[0] = y;
		ypos[1] = y;
		ypos[2] = y;
	}
	return;
}
// --------------------------------------------------	
void plant_1DSOS_RDM::adsorption(){
  //	add_remove_atom((int)(RandN1()*Lx),(int)(RandN1()*Ly),1);
	add_remove_atom((int)(gsl_rng_uniform(gsl_r)*Lx),0,1);
	return;
}
// --------------------------------------------------	
void plant_1DSOS_RDM::migration()
{
	int x,y;
	int xpos[5],ypos[5];
	int i,n;
	int st[4][2];

	migr_site(&x,&y);
	if(lattice_stat[y][x]==4)
		return;

	GetNeighbors(x,y,xpos,ypos);
	
	// Count the number of neighbors with lower height
	for(i=1,n=0;i<N_CL;i++){
		if((lattice[y][x]-lattice[ypos[i]][xpos[i]])>0){
			st[n][0]=ypos[i];
			st[n][1]=xpos[i];
			n++;
		}
	}

	if(n==0) return;
	
	//	n = (int)(RandN1()*n);
	n = (int)(gsl_rng_uniform(gsl_r)*n);
	add_remove_atom(x,y,-1);
	add_remove_atom(st[n][1],st[n][0],+1);

	return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::desorption()
{
	int x,y;
	migr_site(&x,&y);
	add_remove_atom(x,y,-1);
	return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::migr_site(int *migr_x, int *migr_y)
{
	int y,x;
	double rdm;

	//	rdm = RandN1()*migr_ttl;
	rdm = gsl_rng_uniform(gsl_r)*migr_ttl;
	for(y=0;rdm>=0;y++)
		rdm-=p_migr[y];
	y--;
	rdm+=p_migr[y];
	for(x=0;rdm>=0;x++)
		rdm -= MIGR_RT[lattice_stat[y][x]];
	x--;

	*migr_x = x;
	*migr_y = y;
	
	/* FOR DEBUG
	int event[5];
	for(i=0;i<5;i++)
	{
		event[i] = 0;
	}
	if(lattice_stat[y][x]==0)
		event[0]++;
	else if(lattice_stat[y][x] == 1)
		event[1]++;
	else if(lattice_stat[y][x] == 2)
		event[2]++;
	else if(lattice_stat[y][x] == 3)
		event[3]++;
	else if(lattice_stat[y][x] == 4)
		event[4]++;
	*/
	return;
}
// --------------------------------------------------
void plant_1DSOS_RDM::desorp_site(int *desorp_x, int *desorp_y)
{
	int y,x;
	double rdm;

	//	rdm = RandN1()*migr_ttl;
	rdm = gsl_rng_uniform(gsl_r)*desorp_ttl;
	for(y=0;rdm>=0;y++)
		rdm-=p_desorp[y];
	y--;
	rdm+=p_migr[y];
	for(x=0;rdm>=0;x++)
		rdm -= DESORP_RT[lattice_stat[y][x]];
	x--;

	*desorp_x = x;
	*desorp_y = y;

	return;
}
// --------------------------------------------------
double plant_1DSOS_RDM::OneMCTime(double r_d){
	double rd;

	//	rd=RandRand();
	rd = gsl_rng_uniform(gsl_r);
	while ((rd <= 0) || (rd >= 1)){
	  //		rd=RandRand();
	  rd = gsl_rng_uniform(gsl_r);
	}

	return -log(rd)/r_d;
}
// -------------------------------------------------
void plant_1DSOS_RDM::getXCov(const string filename){
	int i;
	int m;
	double h_avg = getH();
	double *xcov_val = new double[Lx];
	for(m=0;m<Lx-1;m++){
		xcov_val[m]=0.0;
		for(i=0;i<Lx-1-m;i++){
			xcov_val[m] += (lattice[0][i]-h_avg)*(lattice[0][i+m]-h_avg);
		}
	}
	ofstream flog(filename.c_str());
	for(m=0;m<Lx-1;m++){
		flog << xcov_val[m];
		if(m!=Lx-1){flog << "\n";}
	}
	return;
}


